{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from [this Jupyter notebook](https://nbviewer.jupyter.org/github/krischer/seismo_live/blob/master/notebooks/Computational%20Seismology/The Finite-Difference Method/fd_ac1d.ipynb) by Heiner Igel ([@heinerigel](https://github.com/heinerigel)), Lion Krischer ([@krischer](https://github.com/krischer)) and Taufiqurrahman ([@git-taufiqurrahman](https://github.com/git-taufiqurrahman)) which is a supplemenatry material to the book [Computational Seismology: A Practical Introduction](http://www.computational-seismology.org/),  additional modifications by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From 1D to 2D acoustic finite difference modelling\n",
    "\n",
    "The 1D acoustic wave equation is very useful to introduce the general concept and problems related to FD modelling. However, for realistic modelling and seismic imaging/inversion applications we have to solve at least the 2D acoustic wave equation.\n",
    "\n",
    "In the class we implement a 2D acoustic FD modelling code based on the 1D code. I strongly recommend that you do this step by yourself starting from this [1D code](https://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/05_2D_acoustic_FD_modelling/lecture_notebooks/1_From_1D_to_2D_acoustic_FD_modelling.ipynb). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite difference solution of 2D acoustic wave equation\n",
    "\n",
    "As derived in [this](https://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/3_Acoustic_medium.ipynb) and [this lecture](https://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/4_2D_1D_elastic_acoustic_approx.ipynb), the acoustic wave equation in 2D with constant density is\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,z,t)}{\\partial t^2} \\ = \\ vp(x,z)^2 \\biggl(\\frac{\\partial^2 p(x,z,t)}{\\partial x^2}+\\frac{\\partial^2 p(x,z,t)}{\\partial z^2}\\biggr) + s(x,z,t) \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "with pressure $p$, acoustic velocity $vp$ and source term $s$. Both second derivatives can be approximated by a 3-point difference formula. For example for the time derivative, we get:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,z,t)}{\\partial t^2} \\ \\approx \\ \\frac{p(x,z,t+dt) - 2 p(x,z,t) + p(x,z,t-dt)}{dt^2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "and equivalently for the spatial derivatives: \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,z,t)}{\\partial x^2} \\ \\approx \\ \\frac{p(x+dx,z,t) - 2 p(x,z,t) + p(x-dx,z,t)}{dx^2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,z,t)}{\\partial x^2} \\ \\approx \\ \\frac{p(x,z+dz,t) - 2 p(x,z,t) + p(x,z-dz,t)}{dz^2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "Injecting these approximations into the wave equation allows us to formulate the pressure p(x) for the time step $t+dt$ (the future) as a function of the pressure at time $t$ (now) and $t-dt$ (the past). This is called an **explicit scheme** allowing the $extrapolation$  of the space-dependent field into the future only looking at the nearest neighbourhood.\n",
    "\n",
    "In the next step, we discretize the P-wave velocity and pressure wavefield at the discrete spatial grid points \n",
    "\n",
    "\\begin{align}\n",
    "x &= i*dx\\nonumber\\\\\n",
    "z &= j*dz\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "with $i = 0, 1, 2, ..., nx$, $j = 0, 1, 2, ..., nz$ on a 2D Cartesian grid.\n",
    "\n",
    "<img src=\"../images/2D-grid_cart_ac.png\" width=\"75%\">\n",
    "\n",
    "Using the discrete time steps\n",
    "\n",
    "\\begin{align}\n",
    "t &= n*dt\\nonumber\n",
    "\\end{align}\n",
    "\n",
    "with $n = 0, 1, 2, ..., nt$ and time step $dt$, we can replace the time-dependent part (upper index time, lower indices space) by\n",
    "\n",
    "\\begin{equation}\n",
    " \\frac{p_{i,j}^{n+1} - 2 p_{i,j}^n + p_{i,j}^{n-1}}{\\mathrm{d}t^2} \\ = \\ vp_{i,j}^2 \\biggl( \\frac{\\partial^2 p}{\\partial x^2} + \\frac{\\partial^2 p}{\\partial z^2}\\biggr) \\ + s_{i}^n. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Solving for $p_{i,j}^{n+1}$ leads to the extrapolation scheme:\n",
    "\n",
    "\\begin{equation}\n",
    "p_{i,j}^{n+1} \\ = \\ vp_{i,j}^2 \\mathrm{d}t^2 \\left( \\frac{\\partial^2 p}{\\partial x^2} + \\frac{\\partial^2 p}{\\partial z^2} \\right) + 2p_{i,j}^n - p_{i,j}^{n-1} + \\mathrm{d}t^2 s_{i,j}^n.\n",
    "\\end{equation}\n",
    "\n",
    "The spatial derivatives are determined by \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,z,t)}{\\partial x^2} \\ \\approx \\ \\frac{p_{i+1,j}^{n} - 2 p_{i,j}^n + p_{i-1,j}^{n}}{\\mathrm{d}x^2} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "and\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,z,t)}{\\partial z^2} \\ \\approx \\ \\frac{p_{i,j+1}^{n} - 2 p_{i,j}^n + p_{i,j-1}^{n}}{\\mathrm{d}z^2}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Eq. (1) is the essential core of the 2D FD modelling code. Because we derived analytical solutions for wave propagation in a homogeneous medium, we should test our first code implementation for a similar medium, by setting\n",
    "\n",
    "\\begin{equation}\n",
    "vp_{i,j} = vp0\\notag\n",
    "\\end{equation}\n",
    "\n",
    "at each spatial grid point $i = 0, 1, 2, ..., nx$; $j = 0, 1, 2, ..., nz$, in order to compare the numerical with the analytical solution. For a complete description of the problem we also have to define initial and boundary conditions. The **initial condition** is \n",
    "\n",
    "\\begin{equation}\n",
    "p_{i,j}^0 = 0, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "so the modelling starts with zero pressure amplitude at each spatial grid point $i, j$. As **boundary conditions**, we assume \n",
    "\n",
    "\\begin{align}\n",
    "p_{0,j}^n = 0, \\nonumber\\\\\n",
    "p_{nx,j}^n = 0, \\nonumber\\\\\n",
    "p_{i,0}^n = 0, \\nonumber\\\\\n",
    "p_{i,nz}^n = 0, \\nonumber\n",
    "\\end{align}\n",
    "\n",
    "for all time steps n. This **Dirichlet boundary condition**, leads to artifical boundary reflections which would obviously not describe a homogeneous medium. For now, we simply extend the model, so that boundary reflections are not recorded at the receiver positions.\n",
    "\n",
    "Let's implement it ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 500 # maximum spatial extension of the 1D model (m)\n",
    "dx   = 1.0 # grid point distance in x-direction\n",
    "\n",
    "tmax = 0.502   # maximum recording time of the seismogram (s)\n",
    "dt   = 0.0010  # time step\n",
    "\n",
    "vp0  = 580.   # P-wave speed in medium (m/s)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 330.0 # receiver position (m)\n",
    "xsrc = 250.0 # source position (m)\n",
    "\n",
    "f0   = 40. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of 2D finite difference with analytical solution\n",
    "\n",
    "In the function below we solve the homogeneous 2D acoustic wave equation by the 3-point spatial/temporal difference operator and compare the numerical results with the analytical solution: \n",
    "\n",
    "\\begin{equation}\n",
    "G_{analy}(x,z,t) = G_{2D} * S \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "with the 2D Green's function:\n",
    "\n",
    "\\begin{equation}\n",
    "G_{2D}(x,z,t) = \\dfrac{1}{2\\pi V_{p0}^2}\\dfrac{H\\biggl((t-t_s)-\\dfrac{|r|}{V_{p0}}\\biggr)}{\\sqrt{(t-t_s)^2-\\dfrac{r^2}{V_{p0}^2}}}, \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "where $H$ denotes the Heaviside function, $r = \\sqrt{(x-x_s)^2+(z-z_s)^2}$ the source-receiver distance (offset) and $S$ the source wavelet.\n",
    "\n",
    "To play a little bit more with the modelling parameters, I restricted the input parameters to dt and dx. The number of spatial grid points and time steps, as well as the discrete source and receiver positions are estimated within this function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": [
     42
    ]
   },
   "outputs": [],
   "source": [
    "# 1D Wave Propagation (Finite Difference Solution) \n",
    "# ------------------------------------------------\n",
    "def FD_1D_acoustic(dt,dx):\n",
    "        \n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction\n",
    "    print('nx = ',nx)\n",
    "        \n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    print('nt = ',nt)\n",
    "    \n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G    = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinate in x-direction\n",
    "\n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - np.abs(x[ir] - x[isrc]) / vp0) >= 0:\n",
    "            G[it] = 1. / (2 * vp0)\n",
    "    Gc   = np.convolve(G, src * dt)\n",
    "    Gc   = Gc[0:nt]\n",
    "    lim  = Gc.max() # get limit value from the maximum amplitude\n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    p    = np.zeros(nx) # p at time n (now)\n",
    "    pold = np.zeros(nx) # p at time n-1 (past)\n",
    "    pnew = np.zeros(nx) # p at time n+1 (present)\n",
    "    d2px = np.zeros(nx) # 2nd space derivative of p\n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------\n",
    "    vp    = np.zeros(nx)\n",
    "    vp    = vp + vp0       # initialize wave velocity in model\n",
    "\n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt) \n",
    "    \n",
    "    # Calculate Partial Derivatives\n",
    "    # -----------------------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # FD approximation of spatial derivative by 3 point operator\n",
    "        for i in range(1, nx - 1):\n",
    "            d2px[i] = (p[i + 1] - 2 * p[i] + p[i - 1]) / dx ** 2\n",
    "\n",
    "        # Time Extrapolation\n",
    "        # ------------------\n",
    "        pnew = 2 * p - pold + vp ** 2 * dt ** 2 * d2px\n",
    "\n",
    "        # Add Source Term at isrc\n",
    "        # -----------------------\n",
    "        # Absolute pressure w.r.t analytical solution\n",
    "        pnew[isrc] = pnew[isrc] + src[it] / dx * dt ** 2\n",
    "                \n",
    "        # Remap Time Levels\n",
    "        # -----------------\n",
    "        pold, p = p, pnew\n",
    "    \n",
    "        # Output of Seismogram\n",
    "        # -----------------\n",
    "        seis[it] = p[ir]   \n",
    "        \n",
    "    # Compare FD Seismogram with analytical solution\n",
    "    # ---------------------------------------------- \n",
    "    # Define figure size\n",
    "    rcParams['figure.figsize'] = 12, 5\n",
    "    plt.plot(time, seis, 'b-',lw=3,label=\"FD solution\") # plot FD seismogram\n",
    "    Analy_seis = plt.plot(time,Gc,'r--',lw=3,label=\"Analytical solution\") # plot analytical solution\n",
    "    plt.xlim(time[0], time[-1])\n",
    "    plt.ylim(-lim, lim)\n",
    "    plt.title('Seismogram')\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.show()            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nx =  500\n",
      "nt =  502\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvEAAAFNCAYAAACJ2DByAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8VFX6x/HPk0JvS5EiCqhIJ4AINopgoUpVEUVhFXZVdHdVFHf9qYvddVdlRV276yKIiEgVaVFARYo0aQKCBFCRJgESUs7vjxmSSTJJBsjkZpLv+/WaF/eee865z00Ok2funHuvOecQEREREZHIEeV1ACIiIiIicnKUxIuIiIiIRBgl8SIiIiIiEUZJvIiIiIhIhFESLyIiIiISYZTEi4iIiIhEGCXxIiIllJklmtk5XschIiInT0m8iEiEM7PLzOxLMztkZvvNbImZXZhfO+dcBefctsKIUUREClaM1wGIiMipM7NKwAzgdmASUAroACR7GVdBMLNo51ya13GIiBRFOhMvIhLZzgdwzk1wzqU554455z5zzq0BMLPfm9kGMztgZnPMrN6JhmbmzOw8/3IPM1tvZofNbJeZ3ecv72xmCWZ2v5n9YmZ7zKyvv/5m/5n/vwb0WdrMXjCz3f7XC2ZWOmD7/f4+dpvZbdlieMfMXjGzWWZ2BLjczHqa2bdm9puZ7TSzRwP6qu9vP8y/7YCZ/dHMLjSzNWZ20MxeCu+PX0TEG0riRUQi22YgzczeNbPuZva7ExvMrC/wV6A/UANYBEzIpZ83gT845yoCzYEFAdtqAWWAM4GHgdeBm4AL8J31fzhgbv3fgIuAVkAc0A54yB9PN+Ae4ArgPKBTkDgGA08AFYHFwBHgZqAK0BO43X9cgdoDDYHrgRf8MVwBNAOuM7Ng+xERiWhK4kVEIphz7jfgMsDhS673mtk0M6sJ/AF4yjm3wTmXCjwJtAo8Gx8gBWhqZpWccweccyuzbXvCOZcCTASqAy865w47574DvgNa+uveCIxxzv3inNsL/B0Y4t92HfC2c+4759xR/7bsPnHOLXHOpTvnkpxz8c65tf71Nfg+hGRPyh/z1/0MX9I/wb//Xfg+uLQO7acpIhI5lMSLiEQ4f5I+1DlXF99Z9Dr4zkjXA170Tys5COwHDN8Z9ewGAD2AHWb2uZldHLBtX8Dc9GP+f38O2H4MqOBfrgPsCNi2w192YtvOgG2By0HLzKy9mS00s71mdgj4I74PEYGyx5JbbCIixYaSeBGRYsQ5txF4B18yvxPfFJkqAa+yzrkvg7Rb5pzrA5wBTMV3keyp2I3vw8MJZ/vLAPYAdQO2nRXsELKtvw9MA85yzlUGXsX3QUREpERTEi8iEsHMrLGZ3Wtmdf3rZwE3AF/jS3gfNLNm/m2VzezaIH2UMrMbzayyf8rMb8Cp3hVmAvCQmdUws+r45tD/z79tEjDMzJqYWTn/tvxUBPY755LMrB2+OfMiIiWekngRkch2GN+FnUv9d3T5GlgH3Ouc+xh4BphoZr/5y7vn0s8QYLu/3h/xXbh6Kh4HlgNrgLXASn8ZzrnZwFhgIbAF+MrfJq/bYd4BjDGzw/iS/lP9hkBEpFgx57J/cykiIhJ+ZtYE3weL0v4Lb0VEJEQ6Ey8iIoXGzPr5p+/8Dt+3BNOVwIuInDwl8SIiUpj+AOwFtuKbd3+7t+GIiEQmTacREREREYkwOhMvIiIiIhJhlMSLiIiIiESYGK8DKMqqVKnizjvvPK/DkCLmyJEjlC9f3uswpAiJtDFx7Bj8vOkQZ6dtJSrHs5Uy7atyHtXOrVyIkRUfkTYmpHBoXEh2gWNixYoVvzrnaoTaVkl8HmrWrMny5cu9DkOKmPj4eDp37ux1GFKERNKY2LYNxlzwCW+kDSDGn8CnWTQ/xl3DT2UbUHr5YtqkfMMb3Mrwg29ww1Xw5JMeBx2BImlMSOHRuJDsAseEme04mbZK4kVESoijR+HeK9fw3sEbifE/kPVorQaUmzOVBi1b0gBITISnOk3miZXdAHjqKWjSBIYM8TBwERHJQXPiRURKiL/8BSpsW00U6QAcq3MO5b79Elq2zKhToQKMWjqQjt0rZJTdcQds317Y0YqISF6UxIuIlADz5sFrr8H/GEJL1rCzRQ/Kzp0OtWrlqBsTAxMnQqNGvvWjiWlMu/LfuPkLCjlqERHJjabTiIgUc8ePw8iRmeutB55H3UkzwXJvU6kSvPsuDLt4I++5G7lgy0p+u6khlX5YA2XKhD9okWIuJSWFhIQEkpKSvA5FClmZMmWoW7cusbGxp9WPp0m8mXUDXgSigTecc09n214a+C9wAbAPuN45t92/7UHgVnxP/LvbOTcnrz7NrCvwD3zfPiQCQ51zW8J9jCIiXnvpJdi0ybdcsSKMHQuWRwJ/Qvv2cP0fqnDuq1sBqPTT96Q89yKxDz0QxmhFSoaEhAQqVqxI/fr1sVD+Q0qx4Jxj3759JCQk0KBBg9Pqy7PpNGYWDYwDugNNgRvMrGm2arcCB5xz5wHPA8/42zYFBgHNgG7Ay2YWnU+frwA3OudaAe8DD4Xz+EREioLDh2HmYys5ix8B+PvfoXbt0Nv/+elaPFX+iYz11CefgYMHCzpMkRInKSmJatWqKYEvYcyMatWqFcg3MF7OiW8HbHHObXPOHQcmAn2y1ekDvOtfngx0Nd9o7wNMdM4lO+d+ALb4+8urTwdU8i9XBnaH6bhERIqMsS86njr4RzZzPuMqPcidQ347qfaVK0ODp0awhXMBKHvsAClPPReOUEVKHCXwJVNB/d69TOLPBHYGrCf4y4LWcc6lAoeAanm0zavP24BZZpYADAGyTN0RESlujhyBlc/OpR3LKEMyI44+T6mUIyfdz7ARsbxYdUzGevrYf8NvJ/dhQESKnujoaFq1apXx2r59O/Hx8VSuXJnWrVvTqFEjOnbsyIwZM057X9u3b6d58+b51nsy24MpLrnkktPed3Hl5Zz4YB9Dsj86MLc6uZUH+1Byos+/AD2cc0vNbBTwL3yJfdYdmo0ARgDUqFGD+Pj4oMFLyZWYmKhxIVkU1THx8cd1uPtw5lSYPT17sHXTpswJ8ich5qY2bBzbiMZsonTSb3x/3wPsGnx9QYZbrBTVMSHeChwXlStX5vDhw57GU7ZsWRYtWpSlbP369Vx88cV8+OGHAKxZs4bBgwfjnDutB1UlJiaSnp6e7zE/+eST3HXXXRnrc+bM8fznFA5JSUnEx8ef3nuFc86TF3AxMCdg/UHgwWx15gAX+5djgF/xJfBZ6p6ol1ufQA1ga0D52cD6/GI8//zznUh2Cxcu9DoEKWKK4phITXWu95krnAPnwKVGxTi3Y8cp93fsmHN/qfh6Rn9HflfHueTkAoy4eCmKY0K8Fzgu1q9f710gfuXLl89RtnDhQtezZ88sZW+++abr27dvjrrx8fEuLi7OxcXFuVatWrnffvvNpaenu/vuu881a9bMNW/e3E2cONE559wPP/zgmjVr5pxz7u2333Z33nlnRj89e/Z0CxcudA888ICLiopycXFxbvDgwVlizK3fhQsXuk6dOrkBAwa4Ro0aucGDB7v09PQC+OmE14nff+CYAJa7k8ilvZxOswxoaGYNzKwUvgtVp2WrMw24xb88EFjgP8hpwCAzK21mDYCGwDd59HkAqGxm5/v7uhLYEMZjExHx1KxZ0H/X2Iz19IHXwdlnn3J/ZcpA1btvYg+++8qXO7Abpkw57ThFxDvHjh3LmErTr1+/XOu1adOGjRs35ih/7rnnGDduHKtWrWLRokWULVuWKVOmsGrVKlavXs28efMYNWoUe/bsCSmep59+mrJly7Jq1SrGjx+fZVte/X777be88MILrF+/nm3btrFkyZKT+ClELs+SeOeb4z4S31n0DcAk59x3ZjbGzK7xV3sTqGZmW4B7gNH+tt8Bk4D1wKfAnc65tNz69JcPBz4ys9X45sSPKqxjFREpbB+O+4UbmJCxHnvvn067z9tGluG1qNsz1hOfHXfafYqI75av4Xrl5UTCvGrVKj7++ONc6/nOn+Z06aWXcs899zB27FgOHjxITEwMixcv5oYbbiA6OpqaNWvSqVMnli1bdjo/HoA8+23Xrh1169YlKioqY25/SeDpfeKdc7OAWdnKHg5YTgKuzaXtE8ATQcpz9Okv/xjIfYSKiBQTu3dD7c/epTTHATgW156y7dqddr+1asEvfYaT8vFjpBDL2uRGXJySAqf5wBIRKdq+/fZbmjRpkqN89OjR9OzZk1mzZnHRRRcxb968XBP+QDExMaSnp2esh3K7xbz6LV26dMZydHQ0qamp+fZXHHg5nUZERMLgnbcdv3dvZKyX/fMfC6zv/nfWph8fcya76PXTGySnK4EXKc7WrFnDY489xp133plj29atW2nRogUPPPAAbdu2ZePGjXTs2JEPPviAtLQ09u7dyxdffEG7bCcR6tevz6pVq0hPT2fnzp188803GdtiY2NJSUnJsa9Q+i1pPD0TLyIiBSs9Hda89AV/ZTMAx8tWotS1Qb/QPCWXXw631uvFwR3Afpg+HQYOLLDuRUqkEE5eF6pFixbRunVrjh49yhlnnMHYsWPp2rVrjnovvPACCxcuJDo6mqZNm9K9e3dKlSrFV199RVxcHGbGs88+S61atbJMcbn00ktp0KABLVq0oHnz5rRp0yZj24gRI2jZsiVt2rTJMi++X79+QfsNNle/pLBQvvYoqRo1auQ2ncKt2KR4i4+PP63bbEnxU5TGxPz5sO2K4QzHdyY+dcTtxPzn5QLdxyOPwBj/beN79ICZMwu0+2KhKI0JKToCx8WGDRuCTlGRkuHE7z9wTJjZCudc21D70HQaEZFi5PXX4c+8wA28z7r6PYkZ/vsC38fQoZnLn34Ku3aUjPmnIiJFiZJ4EZFi4rffYOpUOEp5JnID6Z/MgLYhn9QJWYMGvmk1l7GIV9OHU6X5mXDgQIHvR0REcqckXkSkmJg2DZKTfctxcdCyZfj2NWwY/It7GM4blE/8BSZNCt/OREQkByXxIiLFxMSJmcuDBoV3X/37w8TYmzPWj7763/DuUEREslASLyJSDOzfD0c//YIrmEsUaVx3XXj3V748HOw+iBT/Tc7KrfoStm0L705FRCSDkngRkWLg44/hobRHmctV/Bxbl3N+jA/7PrsPqcGndMss+OijsO9TRER8lMSLiBQDc9/dTWfiAaiW+jOcf37Y99mjB0wvlXmT+KPvTQ77PkWkYH388ceY2Wnfb33o0KFMnpz3e8CTTz6ZZf2SSy45pX09+uijPPfcc6fU9oT4+Hh69eqVZ52DBw/y8suZt+jdvXs3A4vQgzGUxIuIRLhffoHqiz8mCt9zP5IvuRzq1An7fsuVg5Tu12ROqVn7Dfz4Y9j3KyIFZ8KECVx22WVMDLyoJkyyJ/Fffvll2Pd5OrIn8XXq1Mn3g0phUhIvIhLhPvoI+rkpGetlBg8otH33vOl3zCfgSY5TpuReWUSKlMTERJYsWcKbb76ZJYk/8QCigQMH0rhxY2688UZOPBx0zJgxXHjhhTRv3pwRI0aQ/aGh8+fPp1+/fhnrc+fOpX///owePZpjx47RqlUrbrzxRgAqVKiQUe/ZZ5+lRYsWxMXFMXr0aABef/11LrzwQuLi4hgwYABHjx7N83g+/PBDmjdvTlxcHB07dgQgKSmJYcOG0aJFC1q3bs3ChQtztMt+Zr958+Zs376d0aNHs3XrVlq1asWoUaPYvn07zZs3z7Pfd955h/79+9OtWzcaNmzI/fffn89v4dQpiRcRiXALPtxHJz7PLOjbt9D27ZtSk/mh4ch7mhcvEimmTp1Kt27dOP/886latSorV67M2Pbtt9/ywgsvsH79erZt28aSJUsAGDlyJMuWLWPdunUcO3aMGTNmZOmzS5cubNiwgb179wLw9ttvM2zYMJ5++mnKli3LqlWrGD9+fJY2s2fPZurUqSxdupTVq1dnJL79+/dn2bJlrF69miZNmvDmm2/meTxjxoxhzpw5rF69mmnTpgEwbtw4ANauXcuECRO45ZZbSEpKCunn8/TTT3PuueeyatUq/vGPf2TZlle/q1at4oMPPmDt2rV88MEH7Ny5M6T9nSwl8SIiESwxESp9MYMY0gBIbn1RoUylOaFcOUjq1pc0/5+TciuXwJ49hbZ/kWLj0UfBLLTXiBE5248YkbXOo4/mu8sJEyYwyH8/2kGDBjFhwoSMbe3ataNu3bpERUXRqlUrtm/fDsDChQtp3749LVq0YMGCBXz33XdZ+jQzhgwZwv/+9z8OHjzIV199Rffu3fOMY968eQwbNoxy5coBULVqVQDWrVtHhw4daNGiBePHj8+xr+wuvfRShg4dyuuvv05amu89cfHixQwZMgSAxo0bU69ePTZv3pzvzyY/efXbtWtXKleuTJkyZWjatCk7duw47f0FExOWXkVEpFDMnw/XpGVOYSk9qF8etcPjqsE1+HxaJzoTz+ryl9L6l1+gdu1Cj0NEQrdv3z4WLFjAunXrMDPS0tIwM5599lkASpcunVE3Ojqa1NRUkpKSuOOOO1i+fDlnnXUWjz76aNCz2sOGDaN3796UKVOGa6+9lpiYvNNN5xxmlqN86NChTJ06lbi4ON555x3i4+Pz7OfVV19l6dKlzJw5k1atWrFq1aoc032CiYmJIT09PWM9lDP1efUb7GcXDjoTLyISwT77+AhX8VlmQb/CT+K7dYN7Y8ZyJrtoc2QRCdXiCj0GETk5kydP5uabb2bHjh1s376dnTt30qBBAxYvXpxrmxPJbfXq1UlMTMz1Is86depQp04dHn/8cYYOHZpRHhsbS0pKSo76V111FW+99VbGnPf9+/cDcPjwYWrXrk1KSkqOKTjBbN26lfbt2zNmzBiqV6/Ozp076dixY0bbzZs38+OPP9KoUaMs7erXr58xlWjlypX88MMPAFSsWJHDhw8H3Vco/YabkngRkQjlHCR/8ill8f1hPXJOc2jYsNDjqFwZalzenJ/wnX2fPr3QQxCJfI8+6vtPHcrrtddytn/ttax18plOM2HChCwXoAIMGDCA999/P9c2VapUYfjw4bRo0YK+ffty4YUX5lr3xhtv5KyzzqJp06YZZSNGjKBly5YZF7ae0K1bN6655hratm1Lq1atMi4yfeyxx2jfvj1XXnkljRs3zvN4AEaNGkWLFi1o3rw5HTt2JC4ujjvuuIO0tDRatGjB9ddfzzvvvJPlTPmJ496/fz+tWrXilVde4Xz/LXqrVavGpZdeSvPmzRk1alSWNqH0G24WytcMJVWjRo3cpk2bvA5DipgTV+2LnODVmFi5EjZccCM34vujm/63/yPq8TGFHgfASy/BXXf5lrt1g9mzPQmjyND7hAQTOC42bNhAkyZNvA0ojEaOHEnr1q259dZbvQ6lSDrx+w8cE2a2wjnXNtQ+dCZeRCRCzZwJD/E49/BPNtW4jKiB/T2LpXfvzOUFCyCXb6BFpAS44IILWLNmDTfddJPXoRRrurBVRCRCzZwJ22nA89zDhS/eQ6NW3sVSrx60inOUWv0NvY7PIOmiRVRcPQ/yuaBNRIqfFStWeB1CiaB3VxGRCPTLL/DNN77l6Gi4+mpv4wG45hoYvnoAddkF64GvvoIOHbwOS0SkWNJ0GhGRCDR7tu/aNYBLLgH/bZU9dU0fYwa9MtbTP9EVriJ50XWJJVNB/d49TeLNrJuZbTKzLWY2Osj20mb2gX/7UjOrH7DtQX/5JjO7Or8+zecJM9tsZhvM7O5wH5+ISLgsnHKAc9gKQK9e+VQuJG3awFdVM4NJmjwjj9oiJVuZMmXYt2+fEvkSxjnHvn37KFOmzGn35dl0GjOLBsYBVwIJwDIzm+acWx9Q7VbggHPuPDMbBDwDXG9mTYFBQDOgDjDPzM73t8mtz6HAWUBj51y6mZ0R/qMUESl4x49Dlc8msZU/spFGVDp2LzDc67Awg4p9u3L0rbKU4xjldmyArVvh3HO9Dk2kyKlbty4JCQns3bvX61CkkJUpU4a6deuedj9ezolvB2xxzm0DMLOJQB98MylP6AM86l+eDLxkvkd69QEmOueSgR/MbIu/P/Lo83ZgsHMuHcA590sYj01EJGwWL4bLk2YB0JhNuCrHPI4o0xW9yzL/ra70xn8WfsYM+NOfvA1KpAiKjY2lQYMGXochEczL6TRnAjsD1hP8ZUHrOOdSgUNAtTza5tXnufjO4i83s9lmVvhPRBERKQBzPkniCuZlrFvPHh5Gk1XXrjA7KnBKjebFi4iEg5dn4i1IWfaJYbnVya082IeSE32WBpKcc23NrD/wFpDjtglmNgIYAVCjRg3i4+ODBi8lV2JiosaFZFHYYyJhwm+Ux/d48n3V67E2IQESEgpt//nZ2uQy+M63HPvl5yyaOZO08uW9DaqQ6X1CgtG4kOxOZ0x4mcQn4JujfkJdYHcudRLMLAaoDOzPp21u5QnAR/7lj4G3gwXlnHsNeA18T2zVE/ckOz2JUbIrzDHx/ffQbm/m9JTKg/oVufG4YhisvK81bfiW6PRUOhw9Cj17eh1WodL7hASjcSHZnc6Y8HI6zTKgoZk1MLNS+C5UnZatzjTgFv/yQGCB813GPQ0Y5L97TQOgIfBNPn1OBbr4lzsBm8N0XCIiYTNzJvRgVsZ6zDVFZyrNCT16wHQyH+GaPk13qRERKWieJfH+Oe4jgTnABmCSc+47MxtjZtf4q70JVPNfuHoPMNrf9jtgEr4LVj8F7nTOpeXWp7+vp4EBZrYWeAq4rTCOU0SkIK2atJmGbAEgpXR56NjR44hyatwYVtTKnBd/ZOlaD6MRESmePH1iq3NuFgScUvKVPRywnARcm0vbJ4AnQunTX34QKFnf54pIsXL4MFRbmvn2ltrpCmJLl/YwouDM4Ky+F3DPq/9kLlfSo29znvE6KBGRYkZPbBURiRBz58LV6ZlJfNmBRfe8RLceUTzPPayjBbNmB7sXgYiInA4l8SIiEWLe1EQ68XlmQffu3gWTjy5doFQp3/K6dbBzZ971RUTk5CiJFxGJAOnp8Pmnx3iJkWykEUfOawkF8MS/cClfHgJvuDB7tmehiIgUS0riRUQiwMqVsH5vDe7jn3SovpEyS7/wOqR8nfii4Cx+JG3cqzBvXt4NREQkZEriRUQiwIyAuzR27w7RVSt7F0yIevSA3/MmP1KP29fcTtrLr3odkohIsaEkXkQkAsycmbncq1fu9YqShg1hZ91LMtbdnM/g+HEPIxIRKT6UxIuIFHF79sCK5ekAREfDVVd5HFCIzKBRn8b8QH0AYo4ehsWLvQ1KRKSYUBIvIlLEzZ4N87iC+XRhbP1/UsUd8DqkkPXoacwMfETHrByP8RARkVOgJF5EpIiLn7KfTnxOFxZy+7ZRkJrqdUgh69wZ5sZmJvHHp87MvbKIiIRMSbyISBGWnAxR8z4jGt90mqQW7aBGDY+jCl3ZsmCXd+YoZQEotXUjbNvmcVQiIpFPSbyISBH2xRfQJTlzCkqZ/j08jObUXNG7LAvoklmgKTUiIqdNSbyISBE2a3oa3cl8UpL16plH7aKpRw+YReaHj7RpmlIjInK6lMSLiBRRzsGPU5ZTg18BSP5dTWjd2uOoTt4558DGcwI+fMQvhKNHvQtIRKQYUBIvIlJEbd4Mcbsyz1rH9O4BUZH5tt26bz3W0QyA6JRkWLDA44hERCJbZP41EBEpAWbMgB5kzh+P7h158+FP6NED3mYY/2YkQ2vOxnXp6nVIIiIRLcbrAEREJLivpuzhXlYAkB4dQ9SVV3oc0anr0AH6VriXxETgZxj9IzRu7HVUIiKRS2fiRUSKoEOHoOrXmWfhU9p3gMqVPYzo9JQqBYGfQWbq2lYRkdOiJF5EpAj67DOokH6IA1QBoPSAXh5HdPp66sGtIiIFRkm8iEgRNHMmPM89nMEvvHVzPAwa5HVIp61798zlFZ8ncuTTRd4FIyIS4ZTEi4gUMenpmWeqU4ml+Z2doE4db4MqAHXqQNtWqcyiOz+nVaNsz8vhwAGvwxIRiUhK4kVEiphly2DvXt/yGWdA27bexlOQru4ZQzX2UZrjRKWn+eYNiYjISVMSLyJSxMyYkbncI3JvDR9Uz54wk8zJ8U6T40VETkkx+tMgIlI8rJ60ife4ieuZSN9OxWu6Sbt2sKRyZhKfOn22b/6QiIicFE+TeDPrZmabzGyLmY0Osr20mX3g377UzOoHbHvQX77JzK4+iT7/bWaJ4TomEZHTsX07nL95OjcxnoncQI+pw70OqUBFR0Ptnm34iZoAxB7Y65s/JCIiJ8WzJN7MooFxQHegKXCDmTXNVu1W4IBz7jzgeeAZf9umwCCgGdANeNnMovPr08zagv9+bSIiRdC0adCLzPk0sX165lE7MnXvGcVsAm5Vo5vGi4icNC/PxLcDtjjntjnnjgMTgT7Z6vQB3vUvTwa6mpn5yyc655Kdcz8AW/z95dqnP8H/B3B/mI9LROSULfjoAJexOLOgRw/vggmTq6+G2Zb54SRlmubFi4icLC+T+DOBnQHrCf6yoHWcc6nAIaBaHm3z6nMkMM05t6eA4hcRKVAHD0L5xXOIIQ2A43EXQs2aHkdV8KpVg4MXXkkKMQDErl4Be/TWLCJyMmI83LcFKXMh1smtPNiHEmdmdYBrgc75BmU2AhgBUKNGDeLj4/NrIiVMYmKixoVkUVBjYv78M+iWnjm1ZFfrZuwopmPt7BZns/iby7iceAA2Pv88PxWjbx30PiHBaFxIdqczJrxM4hOAswLW6wK7c6mTYGYxQGVgfz5tg5W3Bs4Dtvhm41DOzLb459pn4Zx7DXgNoFGjRq5z586ncmxSjMXHx6NxIYEKaky8/koqfyZzakmDu+6iQZs2p91vUVSlCvzvzZ4ZSfz5W7bRuBj9v9L7hASjcSHZnc6Y8HI6zTKgoZk1MLNS+C5UnZatzjTgFv/yQGCBc875ywf5717TAGgIfJNbn865mc65Ws65+s65+sDRYAm8iIhXjh+HgzMWU439AKTUqAOtWnkcVfjExcHKWj35jYp8RH82NB3gdUgiIhHFszMj3iMGAAAgAElEQVTxzrlUMxsJzAGigbecc9+Z2RhguXNuGvAm8J6ZbcF3Bn6Qv+13ZjYJWA+kAnc659IAgvVZ2McmInKyvvgCrjw6NWM9ZmDf4vWUp2zMoNmAxlQf9ysplOKu32Cs10GJiEQQL6fT4JybBczKVvZwwHISvrnswdo+ATwRSp9B6lQ4lXhFRMLlk6mO+/g4Y9369/MwmsLRt5/x0rhSAEydCi++6EvuRUQkf8X3NI+ISIRwDlZP2Uod/6U9KRWqQKdOHkcVfh07+ubGA+zcCatWeRuPiEgkURIvIuKx1ath0Z7zqMFehpcbT9Tjj0FsrNdhhV1sLPQMeJbVJ1PS4PBh7wISEYkgSuJFRDw2zX9J/yGqcLTvYKL/NNLbgApR377QktW8xnDufqo2PPaY1yGJiEQEJfEiIh6bPDlz+ZprvIvDC1dfDQ1iEhjOG1RN20vKpCm++UUiIpInJfEiIh7auBHWrvUtlymTdXpJSVCxItC1K79REYDYHVth3TpvgxIRiQBK4kVEPPThh/AAT9OZhfTqnkaFEnjvrF4DyzCTgE8vU6Z4F4yISIRQEi8i4qEl47fzNA+ykC6880UD31OfSpi+feETy7yl5vGJH3kYjYhIZFASLyLikY0bocWmDzPWS1/QHEqV8jAib1SvDkmXd+cYZQAotXEtrF/vcVQiIkWbkngREY98+CFczwcZ6zGDr/cwGm/1HlyRGfTKLPjgg9wri4iIkngREa98/b8ttGUFAGkxpaBPH48j8k7fvvBh1KCM9ZT/TdRdakRE8qAkXkTEAxs3QtzmSRnr7sqrMx9fWgJVqwZJXXpk3qVm22b49luPoxIRKbqUxIuIeGD8eLiOzCQ+5saSO5XmhL43lGUqfTMLJk70LhgRkSJOSbyISCFLT4cv395EK1YDkBZbGnr39jgq7wVOqVnMpfxSO87jiEREii4l8SIihWzJEuiy67+ZBT16QKVK3gVURFStCnbVlZzNDjqwmNeO3Oh1SCIiRZaSeBGRQva//6YzhPcy1qNvGeJhNEXL4Fti2cnZAPz3v7q2VUQkN0riRUQKUVISbJy4irokAJBSuRr07JlPq5KjT5/MLyW+/x6WLvU2HhGRokpJvIhIIZoxA75IbMPZ/MizVZ8i5sH7S+QDnnJTtixce23m+n/fdbBvn3cBiYgUUUriRUQK0Xv+WTS7qMuRkaOxB+73NqAi6OaboRq/8jce577XGpI+aLDXIYmIFDkxXgcgIlJS/PQTzJqVuX7TTd7FUpRddhk0qnuUMQkPE5XucPO3wc6dcNZZXocmIlJk6Ey8iEgheestSE31LV92GTRs6G08RVVUFHQddjbz6QqAOee7ylVERDIoiRcRKQRpaTB33GZe4E+czyb++EevIyrahgyBtxmWsZ7yxju6VY2ISAAl8SIiheCzz6DP7pf5E2PZRGOuX/+I1yEVaQ0bwv6O/ThIZQBit2+BhQs9jkpEpOhQEi8iUgjeGXeEobyTsR7T6VLvgokQt44sy3gyH/iUPvYlD6MRESlaPE3izaybmW0ysy1mNjrI9tJm9oF/+1Izqx+w7UF/+SYzuzq/Ps1svL98nZm9ZWax4T4+ERGAhASoPGsCVTgEwPF658EVV3gcVdHXty9Mqn5nxrpN/wR27PAwIhGRosOzJN7MooFxQHegKXCDmTXNVu1W4IBz7jzgeeAZf9umwCCgGdANeNnMovPpczzQGGgBlAVuC+PhiYhkeON1x+1uXMZ6qbtv9129KXmKjYVOtzdl3okLXNPT4ZVXPI5KRKRo8PKvSDtgi3Num3PuODAR6JOtTh/gXf/yZKCrmZm/fKJzLtk59wOwxd9frn0652Y5P+AboG6Yj09EhKNHYemLX9OaVQCkxpaBoUO9DSqCDB8O4+yujPW0V1+HY8c8jEhEpGjw8j7xZwI7A9YTgPa51XHOpZrZIaCav/zrbG3P9C/n2ad/Gs0Q4E/BgjKzEcAIgBo1ahAfHx/yAUnJkJiYqHEhWeQ1Jj75pA43Hxqbsf5Ll8vZvGZNIUVWPOy7+CJ++LI+DdhO9KH9bHzkEX7q0cPrsPKk9wkJRuNCsjudMeFlEm9ByrLfPyy3OrmVB/tmIXufLwNfOOcWBQvKOfca8BpAo0aNXOfOnYNVkxIsPj4ejQsJlNuYSEuDx27ZwnNMyiir8+Tj1GnTphCji3wP/x1evvIO7mYsr8fczl1Dh9O4aQ2vw8qT3ickGI0Lye50xoSX02kSgMDH79UFdudWx8xigMrA/jza5tmnmT0C1ADuKZAjEBHJw9SpcN2P/yCadABSu1wJSuBPWteusKTVSBrwA4+l/pWxE4p2Ai8iUhi8TOKXAQ3NrIGZlcJ3oeq0bHWmAbf4lwcCC/xz2qcBg/x3r2kANMQ3zz3XPs3sNuBq4AbnXHqYj01ESjjn4M3H92S9reRDD3oXUAQzg7/8tSxp/i+Px42Dw4c9DkpExGOeJfHOuVRgJDAH2ABMcs59Z2ZjzOwaf7U3gWpmtgXf2fPR/rbfAZOA9cCnwJ3OubTc+vT39SpQE/jKzFaZ2cOFcqAiUiLFx8OOVftZxoUAHG/THvQ1+inr3x/OO8+3fOAAvP66t/GIiHjNyznxOOdmAbOylT0csJwEXJtL2yeAJ0Lp01/u6bGKSMnhHPz1r7CeZnRgMc/0WsT9D0b7TinLKYmOhlGj4A9/8K3Pfno1d+//mJjHH/UwKhER7+hGxSIiBWz6dPjaf/+sUqVg0LgOcMkl3gZVDNx8M9Q6I53JDGDu3lbEPPF3WBT0HgUiIsWekngRkQKUlgZ/+1vm+u23w9lnexdPcVKmDDzwYBQHqZJRljpqtO+rDxGREkZJvIhIAZo4Ec5c9ymlSKZ8ed+0Gik4t98O79T5G8mUAiBm6Zfw4YceRyUiUvjyTeLNrJyZ/Z+Zve5fb2hmvcIfmohIZDl2DMbfv5pZ9GAjjXnrygmcUUNniQtS6dLwh2fOYSx3Z5Sl3vcAJCV5GJWISOEL5Uz820AycLF/PQF4PGwRiYhEqCefcPx59yiicDRgO/2Ova+LWcNg8GCY1vxv/Eo1AGJ2boexY/NuJCJSzISSxJ/rnHsWSAFwzh0j+BNTRURKrI0bYfvTE7mKuQCkWxSx/3za46iKp6goeOi5KjzC3zPKUv/+OOzZ42FUIiKFK5Qk/riZlQUcgJmdi+/MvIiI4LuucvRtv/KvtMwpHnbHHdCsmYdRFW9XXQU7u/+B9TQBIOboYdKG/0EXuYpIiRFKEv8IvgcqnWVm44H5wP1hjUpEJILMmVOT/kvuoQa/AnC81lnYU096HFXxZgYvvRrDvWVeziiLnjkd3n/fw6hERApPvkm8c24u0B8YCkwA2jrn4sMblohIZNiyBTa+sI2beS+jrNQbr0DFih5GVTKcfTZ0e7oz47gjoyzlvgchJcXDqERECkeuSbyZtTnxAuoBe4DdwNn+MhGREu34cbh7wC5eSx6WUZZy7Q3Qs6eHUZUsI0fChxc8ww/UZyGduaZSPEdTYr0OS0Qk7GLy2PZP/79lgLbAanwXtLYElgKXhTc0EZGi7aEHUvjbmus4g70ApFSrSexLL3gcVckSHQ3/frsCV7b9gm3Hz8RtjuLOO+Htt72OTEQkvHI9E++cu9w5dzmwA2jjnGvrnLsAaA1sKawARUSKog8+gJ0vTOZSvgT8d6P56AM44wyPIyt5WrSAB146C+f/k/bOO/DWW97GJCISbqFc2NrYObf2xIpzbh3QKnwhiYgUbfHxcPPNMJFBDOc1kq009uST0KmT16GVWLfdBkOGZK7fcQes+dc82LDBu6BERMIolCR+g5m9YWadzayT/8mtelcUkRLpu++gb1/ffHgwljQZzuJxb2H3j/I6tBLNDF555cRdPR13Jz9Ls3uv5ninK2DrVq/DExEpcKEk8cOA74A/AX8G1vvLRERKlDVr4OquqRw65FuvVQtmz4boJnV8TyAST5UvD598Au2rb+NhxhBNOqX27ialYxfYscPr8EREClQot5hMcs4975zr538975xLKozgRESKiiVL4KWLxzP759ZUZR8VKsCsWVCvnteRSaBzz4XX5p/LDeWnc4wyAMTu/pHkSy73PVZXRKSYyDeJN7MfzGxb9ldhBCciUhRM/tAxt/MTvHb0JlqwjhnRfZg1JYnWrb2OTIJp2RJGz7mc60tNJZlSAJTe/QMpbS+CuXM9jk5EpGCE8v1vW+BC/6sDMBb4XziDEhEpCpKSYNSwX7HrBvJo6kMZ5a3qH6RD8wMeRib5ufRSGDXvam6uMIUjlAMg9sgh0q/uTtpTz0JamscRioicnlCm0+wLeO1yzr0AdCmE2EREPLNkCdzbZBb3vNOCAUzJKD96URfKLl8MtWt7GJ2EokMH+Puyngyqs4gEzgQgyqUR/dcHONa2A2za5HGEIiKnLpTpNG0CXm3N7I+AnicuIsVSQgI82Gsthy7rwbjtPanNTxnbkn9/O+U+nw1VqngYoZyMxo3hjZVtuKv9Mr6mfUZ52VVfsa7HKPbt8zA4EZHTEMp0mn8GvJ4C2gDXhTMoEZHCtm4djBqwjSVnD+KJmXH0YHbGtqOVauJmzKT0my9DqVIeRimnomZNmLykNoufXMSjUWM4TiwAt2x7hHPPhYcegp07PQ5SROQkhZLE33ri6a3OuSudcyOA4+EOTEQk3HbvhhdfhPbtfU/9nDClFAPcZKJwAKRjJF73e8ptWYv17OFxtHI6oqPhvgdj6bfy/7gtbjlj+D9WcgGHDsETT0CDBnBzj1/Z1bYPx94YDwd0zYOIFG2hJPGTQyw7aWbWzcw2mdkWMxsdZHtpM/vAv32pmdUP2Pagv3yTmV2dX59m1sDfx/f+PnU6TaQESUuDTRvSmTF2G+/2nMSEaiO59czZ/PnP8M03vjq7qMvH9ANg34XdiFqzmgofvAk1angYuRSkuDh499uWtJgyhvPPzyxPS4MzZr/DmSumUXb4TaRVq8GOBp3ZccvDJH40B827EZGiJia3DWbWGGgGVDaz/gGbKoH/5runwcyigXHAlUACsMzMpjnn1gdUuxU44Jw7z8wGAc8A15tZU2CQP746wDwzO/F2nFufzwDPO+cmmtmr/r5fOd3jEBHvpaTAwYO+k6eHd+wnae33HNy2n+NbdxK9czulf9pB5QM/0DR9HY04nNHuN1L4lO4AxMbCNddAw4FjoMlDVIuL8+pwJMzMoF8/6N0bpk+Hl16CBQsct/FGRp1ol0a97Z/D9s/hv76yg7HV+bVGE46c1Zifuw+jVKeLqV3bd4lEpUpQxh3Dypbx7UBEJMxyTeKBRkAvoArQO6D8MDC8APbdDtjinNsGYGYTgT74ngh7Qh/gUf/yZOAlMzN/+UTnXDLwg5lt8fdHsD7NbAO+O+oM9td5199vnkn88R2H+KJxPofq+9adBU1HklAt6x/965beS9njhwKr5WpG3EPsq1gfF1Bx2OLf57vfEyZc+C+Olc682K5c8gFu+OYv+ew1s6s3Ln0nS1n1xO30W/VISO2PlKrC+AtfzOzPQb3939Jtw7/y3zHwa4X6TI57LMumZj/Np+PWt7LUy832qq2Z1eS+LD+7C3d+RLudH+Woay5nZ+tqdmXhubdlKbt8y+s0/3l+0P2lJyfxZelXM9a/rjuQpWcNzFLnmg3PcM6BFXkH7jfvnBGsrXkFQMYx3LjmAWoeyfo4Bgv4QQQexdRGo9n6u7ZZ6v5xxW1UOB7adIB3mz3LzxXOzVL2wNf9ctTL7dcwrvWbJJaqmrFePnk/d30bbOwG7+HJdp9kWa+VuIVbv7sns5UDl5ZOdFoy0anJxKYlE5uWRGx6Mr9YTa4q/TmJiZntb2Y67zI0l2iz6sQXdOkC118PAwdC1aoATUJqK5EvJsaXzPfrBxs3GvP+M5O5H3xAuz2f0J5vctSvkvIrVXYvgt2LGLr0Et7l4izbN9KahmwmMaoSDa0i22IqkBYVS1pUKVKjS5EeHUua/9//xv2LPZUbY0bG664vbwAD8BW4ExuCrE+6KOt7ftnkgwz8ZlTIx/5eh9ezrFc7vJ0eq54Mqe2xUpWZ3P4fWcrq7lvN5evHhdR+X4V6zGr9tyxlDfd8wUVbQrt7dULVlixsNjJLWdyOabT8cUZI7b+v1YGvGw7JUnbR9+/R8KdFIbVffXZv1tTrnaWsy3f/5sz9a/Ntm3T0KM+2jOL72h2zlPf89nGqJv4Y0v4XNB3Jrmots5QNXDoqI9/Iz8xWf2N/xYCn1DnHkMV/CKktwOR2zwYZe/eH3P69Dq9lWc8x9vL4DHysVGU+CjL2Ooc49vbnMvbahzj2duUy9lqEOPa2BBl7jb5ZQ1Ty++wvf3ZIfWThnMvzBVycX51TeQEDgTcC1ocAL2Wrsw6oG7C+FagOvATcFFD+pr+/oH3622wJKD8LWJdfjE0o55w/h8jv1YtpOYp3Uyvk9m1YnqM41LYOXC12Zymqza6Tap+9qA3LQ267i9o5insxLeT2y2mTo3gEr4bcfhq9chQ/wiMht3+VETmK/8PwkNs/wiM5iqfRK+T2w/lPjuLltAm5fUkeewnUyVGc39g7GFvNfV+vi/v+hofc0elzXUFYuHBhgfQjRcP27c59OHa3+0+Xie79aiPdcrvAHaFslnF0JXNyDK9d1A557F7AsmxF6Z7+v7uAZSG3Dfae35tPQm6/jAtyFJ/Me/4n9M5R/CgPh9z+dN/zH+bRHMXT6Rly+2Dv+StoHXL7nkzPUXwy7/mtWZGtSGMv1PbhHHu+XIjlzoWeS+c1neZ+59yzwGAzuyFI8n/3yX9kyLqLIGUuxDq5lQeb459X/ZxBmY0ARgA08T8gRESKptIkA2DmqFAhlYoVU4mNqcz6X+NIKluJo1XO4FjNWqTVrUG5plWJaVqL49WrZ0x3SACIjz/tOBITE4kvgH6k6KjeAqq3qAkM4GDaQKbsKsWBNb/h1u+i4u4dVCpXixaJB9m3rzRHj0Zz5EgMMSmpIfd/nKyXZVm+39dm5fI6XSkiJYI5F/yNw8x6O+emm9ktwbY75949rR2bXQw86py72r/+oL/fpwLqzPHX+crMYoCfgBrA6MC6J+r5m+XoE3ga2AvUcs6lZt93bhpUOdu9e83/5X8swK5mV5FYrV6W8nO/Hk9MyrF82wPsaN2X5IrVff3535sbLnqTvL5XCpx2ua3dINJKZ37oiE4+SoNlk3LUy83WDkOzrJc+/Ct1V88MstOcRWmxZdnRPutdR8vv30nNjfH57teA4xWqsqtVzyzlFX/6nmpbM7/OzusYjv7uTH5p2jlLvSo/rqHyznXB9xmVtbPDZ5zL/vPaZSmruuUbKu7NOp3lhISEBOrWrZux/tvZzTl0VvMsdapvWETZg3uy7Th4/AcaXMCRWpnTWczgjDVzKXXkYI66LuAHcWJpX6NLSKpaJ0u92itmEJWSnOfXkifsbd6F1Iq/y1JW5+spIc/r/bl1N9JLl81Yj0o+Rq3Vc4LWzZ54mMFP7ftkKYtJPEj19V9k9hcFMbFGTPnSRJcrTUyFMsSUL01sxTLEVixDmcb1qVTJV88r8fHxdO7c2bsApEhITobf9qeSuOcwyxYsokWD80g9lkJaUgqpx46Tfuw4aUkppCcd5+fzO5BStlLm+bh0x9lfTsw4QefSs56wy76+pd1g0kpl/r+LST7Cud+8H3KsmzpknSpa5vBe6n/7cUhtU0uVY8tFN2Upq/Drduqu/yyk9kkVqrO9Tf8sZVV2r6fWlsUhtU+sVo+EZln/fFffvpzqP64Mqf3B2k34qWGHLGW1Nn9BlZ82htT+13oX8Gu9C7KU1V33KRX25z8d5qeff4LLruVg7azT9uqv+IgyR0K7eDqh2dU58o3zvnov5Hzjh9b9M/INAJyj8aLXc2+Qzfftb8qSb8QkH+G8peNDbr+x44gs64FjL5eUNENxHHsp8/9Hi9ijJFWoztX/GbDCOdc2l+Y55JrEh5s/Kd8MdAV2AcuAwc657wLq3Am0cM790X9ha3/n3HVm1gx4H988+DrAfKAhvpQlaJ9m9iHwkcu8sHWNc+7lvGJs1KiR26Qn+kk2StgkO40JyU5jQoLRuJDsAseEmZ1UEp/XdJrp5DLlBMA5d81JxBisfaqZjQTmANHAW/5kewy+OUHT8M11f89/4ep+fHekwV9vEr6LYFOBO51zaf64c/Tp3+UDwEQzexz41t+3iIiIiEjEyevuNM+Fe+fOuVnArGxlDwcsJwHX5tL2CeCJUPr0l28j8w42IiIiIiIRK9ck3jn3+Yll/4ORGuM7M7/JOacntoqIiIiIeCSvM/EAmFlP4FV8t3c0oIGZ/cE5NzvcwYmIiIiISE75JvHAP4HLnXNbAMzsXGAmoCReRERERMQDodyU7ZcTCbzfNuCXMMUjIiIiIiL5COVM/HdmNguYhG9O/LXAMjPrD+CcmxLG+EREREREJJtQkvgywM9AJ//6XqAq0BtfUq8kXkRERESkEOWbxDvnhhVGICIiIiIiEppQ7k7TALgLqB9Y/3Qf9iQiIiIiIqcmlOk0U/E93XQ6kB7ecEREREREJD+hJPFJzrmxYY9ERERERERCEkoS/6KZPQJ8BiSfKHTOrQxbVCIiIiIikqtQkvgWwBCgC5nTaZx/XUREREREClkoSXw/4Bzn3PFwByMiIiIiIvkL5Ymtq4Eq4Q5ERERERERCE8qZ+JrARjNbRuaceOec6xO+sEREREREJDehJPGPBCwbcBlwQ3jCERERERGR/OQ7ncY59zlwCOgJvAN0BV4Nb1giIiIiIpKbXM/Em9n5wCB8Z933AR8A5py7vJBiExERERGRIPKaTrMRWAT0ds5tATCzvxRKVCIiIiIikqu8ptMMAH4CFprZ62bWFd+ceBERERER8VCuSbxz7mPn3PVAYyAe+AtQ08xeMbOrCik+ERERERHJJpQLW48458Y753oBdYFVwOiwRyYiIiIiIkGF8rCnDM65/c65/zjnuoQrIBERERERydtJJfEFxcyqmtlcM/ve/+/vcql3i7/O92Z2S0D5BWa21sy2mNlYM7O8+jWzG81sjf/1pZnFFc6RioiIiIgUPE+SeHzTceY75xoC8wkyPcfMquJ70FR7oB3wSECy/wowAmjof3XLp98fgE7OuZbAY8Br4TgoEREREZHC4FUS3wd417/8LtA3SJ2rgbn+KTwHgLlANzOrDVRyzn3lnHPAfwPaB+3XOfelvw+Ar/HN7RcRERERiUh53Sc+nGo65/YAOOf2mNkZQeqcCewMWE/wl53pX85eHmq/twKzcwvMzEbgO8tPjRo1iI+PD+mApORITEzUuJAsNCYkO40JCUbjQrI7nTERtiTezOYBtYJs+luoXQQpc3mUhxLT5fiS+Mtyq+Ocew3/dJtGjRq5zp07h9K1lCDx8fFoXEggjQnJTmNCgtG4kOxOZ0yELYl3zl2R2zYz+9nMavvPltcGfglSLQHoHLBeF9/96hPIOh2mLrDbv5xrv2bWEngD6O6c23cKhyQiIiIiUiR4NSd+GnDibjO3AJ8EqTMHuMrMfue/oPUqYI5/usxhM7vIf1eamwPaB+3XzM4GpgBDnHObw3FAIiIiIiKFxask/mngSjP7HrjSv46ZtTWzN8B3T3p8d5JZ5n+N8ZcB3I7vrPoWYCuZc9yD9gs8DFQDXjazVWa2PMzHJyIiIiISNp5c2OqfztI1SPly4LaA9beAt3Kp1/wk+r0tsF8RERERkUjm1Zl4ERERERE5RUriRUREREQijJJ4EREREZEIoyReRERERCTCKIkXEREREYkwSuJFRERERCKMkngRERERkQijJF5EREREJMIoiRcRERERiTBK4kVEREREIoySeBERERGRCKMkXkREREQkwiiJFxERERGJMEriRUREREQijJJ4EREREZEIoyReRERERCTCKIkXEREREYkwSuJFRERERCKMkngRERERkQijJF5EREREJMIoiRcRERERiTBK4kVEREREIownSbyZVTWzuWb2vf/f3+VS7xZ/ne/N7JaA8gvMbK2ZbTGzsWZmofRrZheaWZqZDQzvEYqIiIiIhI9XZ+JHA/Odcw2B+f71LMysKvAI0B5oBzwSkJS/AowAGvpf3fLr18yi/7+9+4+1u77rOP58bc38tTLLuEBtl7Gw0oWxwULH0GV6J22pM7aYYZxOKcsIGcLwZwIRExT8o4gyXZxorU2KcWHK5uhcWFfQu+mE0eqA8kPWMheorcBos7UTp5W3f9xvw+nh3HsP6z3n9Nv7fCQn9/v5fD/ncz7f9J3TV7/9nHOBm4Atg7ggSZIkaVhGFeLXAJua403ART3GXAhsrap9VbUf2AqsSrIQOKGq7q2qAm7reP50834Y+CTwzKxeiSRJkjRkowrxp1TVXoDm58k9xiwCnupo7276FjXH3f1TzptkEfDTwJ/O4jVIkiRJIzFvUBMnuRs4tcep6/qdokdfTdM/nT8Erqmq/2u2z0/9osnlTG7VYWxsjImJiZlXqjnl4MGD1oWOYE2omzWhXqwLdTuamhhYiK+q5VOdS/J0koVVtbfZHtNri8tuYLyjvRiYaPoXd/XvaY6nmncZcHsT4E8C3pPkUFV9use61wPrAZYuXVrj4+PdQzTHTUxMYF2okzWhbtaEerEu1O1oamJU22k2A4e/bWYtcGePMVuAlUkWNB9oXQlsabbJHEhyfvOtNJd0PL/nvFX1hqo6rapOA+4AfqlXgJckSZLaYFQhfh2wIslOYEXTJsmyJBsAqmofcCOwrXnc0PQBXAFsAHYBTwB3TTevJEmSdDwZ2Haa6VTVc8AFPfq3A5d1tDcCG6cYd1a/83aNufTlr1iSJEk6dvgbWyVJkqSWMcRLkiRJLWOIlyRJklrGEC9JkiS1jCFekiRJahlDvCRJktQyhnhJkqN8pcYAAAtVSURBVCSpZQzxkiRJUssY4iVJkqSWMcRLkiRJLWOIlyRJklrGEC9JkiS1jCFekiRJahlDvCRJktQyhnhJkiSpZQzxkiRJUssY4iVJkqSWMcRLkiRJLWOIlyRJklrGEC9JkiS1jCFekiRJahlDvCRJktQyhnhJkiSpZUYS4pOcmGRrkp3NzwVTjFvbjNmZZG1H/7lJdiTZleSjSTLTvEnGkzyQ5JEkXxj8VUqSJEmDMao78dcC91TVEuCepn2EJCcC1wPvAM4Dru8I5bcClwNLmseq6eZN8oPAnwCrq+rNwM8M6LokSZKkgRtViF8DbGqONwEX9RhzIbC1qvZV1X5gK7AqyULghKq6t6oKuK3j+VPN+/PAp6rqSYCqema2L0iSJEkallGF+FOqai9A8/PkHmMWAU91tHc3fYua4+7+6eY9A1iQZCLJvyS5ZNauRJIkSRqyeYOaOMndwKk9Tl3X7xQ9+mqa/unMA84FLgC+D7g3yX1V9dWXvGhyOZNbdRgbG2NiYqLP5WquOHjwoHWhI1gT6mZNqBfrQt2OpiYGFuKravlU55I8nWRhVe1ttsf02t6yGxjvaC8GJpr+xV39e5rjqebdDXyjqr4NfDvJF4GzgZeE+KpaD6wHWLp0aY2Pj3cP0Rw3MTGBdaFO1oS6WRPqxbpQt6OpiVFtp9kMHP62mbXAnT3GbAFWJlnQfKB1JbCl2SZzIMn5zbfSXNLx/KnmvRN4V5J5Sb6fyQ/LPjbbFyVJkiQNw6hC/DpgRZKdwIqmTZJlSTYAVNU+4EZgW/O4oekDuALYAOwCngDumm7eqnoM+BzwEHA/sKGqHh70RUqSJEmDMLDtNNOpqueY3J/e3b8duKyjvRHYOMW4s/qdtzl3M3Dzd79qSZIk6djgb2yVJEmSWsYQL0mSJLWMIV6SJElqGUO8JEmS1DKGeEmSJKllDPGSJElSyxjiJUmSpJYxxEuSJEktY4iXJEmSWsYQL0mSJLWMIV6SJElqGUO8JEmS1DKGeEmSJKllDPGSJElSyxjiJUmSpJYxxEuSJEktY4iXJEmSWsYQL0mSJLWMIV6SJElqGUO8JEmS1DKGeEmSJKllDPGSJElSy4wkxCc5McnWJDubnwumGLe2GbMzydqO/nOT7EiyK8lHk2S6eZO8JslnkjyY5JEkHxjOlUqSJEmzb1R34q8F7qmqJcA9TfsISU4ErgfeAZwHXN8R9m8FLgeWNI9VM8x7JfBoVZ0NjAN/kORVA7guSZIkaeBGFeLXAJua403ART3GXAhsrap9VbUf2AqsSrIQOKGq7q2qAm7reP5U8xYwv7lj/2pgH3Bolq9JkiRJGop5I3rdU6pqL0BV7U1yco8xi4CnOtq7m75FzXF3/3Tz/jGwGdgDzAd+tqpemK2LkSRJkoZpYCE+yd3AqT1OXdfvFD36apr+6VwIPAD8OHA6sDXJP1bVt17yosnlTG7VYWxsjImJiT6Xq7ni4MGD1oWOYE2omzWhXqwLdTuamhhYiK+q5VOdS/J0koXN3fKFwDM9hu1mcv/6YYuBiaZ/cVf/nuZ4qnk/AKxrtt/sSvLvwJuA+3usez2wHmDp0qU1Pj7ePURz3MTEBNaFOlkT6mZNqBfrQt2OpiZGtSd+M3D422bWAnf2GLMFWJlkQfOB1pXAlma7zIEk5zd73C/peP5U8z4JXACQ5BRgKfC12b0kSZIkaThGFeLXASuS7ARWNG2SLEuyAaCq9gE3Atuaxw1NH8AVwAZgF/AEcNd08zbz/EiSHUx+a801VfWNwV6iJEmSNBgj+WBrVT1Hc2e8q387cFlHeyOwcYpxZ72MefcweSdfkiRJaj1/Y6skSZLUMoZ4SZIkqWUM8ZIkSVLLGOIlSZKkljHES5IkSS2Tyd9/pF6SHAAeH/U6dMw5CfArStXJmlA3a0K9WBfq1lkTr6+qsX6fOJKvmGyRx6tq2agXoWNLku3WhTpZE+pmTagX60LdjqYm3E4jSZIktYwhXpIkSWoZQ/z01o96ATomWRfqZk2omzWhXqwLdfuua8IPtkqSJEkt4514SZIkqWUM8UCSVUkeT7IrybU9zn9Pkk8057+c5LThr1LD1EdN/GiSf01yKMnFo1ijhq+Puvi1JI8meSjJPUleP4p1anj6qIkPJdmR5IEk/5TkzFGsU8MzU010jLs4SSXx22rmgD7eKy5N8mzzXvFAkstmmnPOh/gkrwQ+BvwEcCbwcz3eZD8I7K+qNwIfAW4a7io1TH3WxJPApcDHh7s6jUqfdfEVYFlVvRW4A/i94a5Sw9RnTXy8qt5SVecwWQ+3DHmZGqI+a4Ik84GrgS8Pd4UahX7rAvhEVZ3TPDbMNO+cD/HAecCuqvpaVf0PcDuwpmvMGmBTc3wHcEGSDHGNGq4Za6Kqvl5VDwEvjGKBGol+6uIfquq/muZ9wOIhr1HD1U9NfKuj+QOAH0Q7vvWTKQBuZPIfdf89zMVpZPqti5fFEA+LgKc62rubvp5jquoQ8E3gtUNZnUahn5rQ3PNy6+KDwF0DXZFGra+aSHJlkieYDG1XD2ltGo0ZayLJ24DXVdXfDXNhGql+//54b7Md844kr5tpUkM89Lqj3n2npJ8xOn74561e+q6LJL8ALANuHuiKNGp91URVfayqTgeuAX5r4KvSKE1bE0leweS23F8f2op0LOjnveIzwGnNdsy7eXEHyJQM8ZP/Gur8185iYM9UY5LMA14D7BvK6jQK/dSE5p6+6iLJcuA6YHVVfWdIa9NovNz3ituBiwa6Io3aTDUxHzgLmEjydeB8YLMfbj3uzfheUVXPdfyd8efAuTNNaoiHbcCSJG9I8irgfcDmrjGbgbXN8cXA35dfsH8866cmNPfMWBfNf5P/GZMB/pkRrFHD1U9NLOlo/iSwc4jr0/BNWxNV9c2qOqmqTquq05j87Mzqqto+muVqSPp5r1jY0VwNPDbTpPNmdYktVFWHklwFbAFeCWysqkeS3ABsr6rNwF8Af5lkF5N34N83uhVr0PqpiSRvB/4WWAD8VJLfqao3j3DZGrA+3ytuBl4N/E3z2fcnq2r1yBatgeqzJq5q/nfmf4H9vHhDSMehPmtCc0yfdXF1ktXAISaz5qUzzetvbJUkSZJaxu00kiRJUssY4iVJkqSWMcRLkiRJLWOIlyRJklrGEC9JkiS1jCFekiRJahlDvCTNcUlem+SB5vGfSf6jo/3PA3rNtyXZMM35sSSfG8RrS9LxYM7/sidJmuuq6jngHIAkvw0crKrfH/DL/ibwu9Os6dkke5O8s6q+NOC1SFLreCdekjSlJAebn+NJvpDkr5N8Ncm6JO9Pcn+SHUlOb8aNJflkkm3N45095pwPvLWqHmzaP9Zx5/8rzXmATwPvH9KlSlKrGOIlSf06G/hl4C3ALwJnVNV5wAbgw82YPwI+UlVvB97bnOu2DHi4o/0bwJVVdQ7wLuD5pn9705YkdXE7jSSpX9uqai9AkieAzzf9O4B3N8fLgTOTHH7OCUnmV9WBjnkWAs92tL8E3JLkr4BPVdXupv8Z4Idm/zIkqf0M8ZKkfn2n4/iFjvYLvPj3ySuAH66q55na88D3Hm5U1boknwXeA9yXZHlV/VszZrp5JGnOcjuNJGk2fR646nAjyTk9xjwGvLFjzOlVtaOqbmJyC82bmlNncOS2G0lSwxAvSZpNVwPLkjyU5FHgQ90Dmrvsr+n4AOuvJHk4yYNM3nm/q+l/N/DZYSxaktomVTXqNUiS5pgkvwocqKrpviv+i8Caqto/vJVJUjt4J16SNAq3cuQe+yMkGQNuMcBLUm/eiZckSZJaxjvxkiRJUssY4iVJkqSWMcRLkiRJLWOIlyRJklrGEC9JkiS1zP8DSYn0uYVExucAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dx   = 1.0     # grid point distance in x-direction (m)\n",
    "dt   = 0.0010  # time step (s)\n",
    "FD_1D_acoustic(dt,dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* How to implement a 2D acoustic FD modelling code based on an existing 1D code"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
